Regional Cumulative Impacts

bc  = readOGR('shapefiles','ohibc_rgn', verbose = F)%>%
      spTransform(CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
ne  = readOGR('shapefiles','ohine_rgn', verbose = F)%>%
      spTransform(CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))

#baffin bay and davis strait
bbds_map <- readOGR(file.path(dir_M,'git-annex/mci-report/shapefiles'),'baffin', verbose = F)

bbds <- spTransform(bbds_map,CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
bbds@data <- bbds@data%>%dplyr::select(-id)

mci_rgns = readOGR(dsn = file.path(dir_M,'git-annex/mci-report/shapefiles'),layer = 'mci_rgns', verbose=F)%>%
    spTransform(CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
mci_rgns@data <- mci_rgns@data%>%dplyr::select(-ohi_rgn)

mci_rgns = rbind(mci_rgns,bbds)

if(!file.exists(file.path(dir_M,"git-annex/mci-report/three_nm.tif"))){

#3nm shapefile to use when clipping coastal pressures
threenm <- readOGR(dsn = file.path(dir_M,'git-annex/globalprep/spatial/d2014/data'),layer = 'regions_offshore3nm_mol')

#remove all polygons outside of our regions
t <- threenm[as.character(threenm@data$rgn_name) %in% c(as.character(mci_rgns@data$rgn_nam),"United States","Russia","Canada","Greenland"),]

#rasterize to use for masking (I think this is faster than masking with shapefile)
three_ras <- rasterize(t,ocean,t@data$sp_id,progress='text', filename = file.path(dir_M,'git-annex/mci-report/three_nm.tif'),overwrite=T)
}


three_ras <- raster(file.path(dir_M,'git-annex/mci-report/three_nm.tif'))
three_ras_bc <- crop(three_ras,bc)%>%mask(bc)
three_ras_ne <- crop(three_ras,ne)%>%mask(ne)
#three_ras_gd <- crop(three_ras,bbds)%>%mask(bbds)

#read in from Mazu location: marine_threats/impact_layers_2013_redo/global_impact_model_2013/normalized_by_two_time_periods/averaged_by_num_ecosystems/by_threat

layers <- list.files(file.path(dir_M,'marine_threats/impact_layers_2013_redo/global_impact_model_2013/normalized_by_two_time_periods/averaged_by_num_ecosystems/by_threat'),full.names=T, pattern = '.tif$')[c(1:5,7,9:12,14,15)]

Regions

Reporting regions

getvals <- function(layer){
  
    if(basename(layer) %in% c("artisanal_fishing_combo.tif","inorganic_combo.tif","plumes_fert_combo.tif","plumes_pest_combo.tif")){
      
       vals = raster(layer)%>%
              crop(three_ras)%>%
              mask(three_ras)%>%
              raster::extract(mci_rgns,method = 'simple',progress='text')%>%
              setNames(mci_rgns@data$rgn_nam)
      
    }else{
      
       vals = raster(layer)%>%
              crop(ocean)%>%
              raster::extract(mci_rgns,method = 'simple',progress='text')%>%
              setNames(mci_rgns@data$rgn_nam)
    }  
return(vals)
}

#check to see if file exists.
if(!file.exists(file.path(dir_M,"git-annex/mci-report/rgn_prs_values.RData"))){

  #list of all values within each region for each pressure
mci_prs_vals <- mclapply(layers,getvals,mc.cores=12)%>%
               setNames(basename(layers))%>%
               setNames(str_replace(names(.),"_combo.tif",""))

  save(mci_prs_vals, file=file.path(dir_M,"git-annex/mci-report/rgn_prs_values.RData"))
}
  #vals is the list, let's try to use map to calculate mean, max
  
  load(file.path(dir_M,"git-annex/mci-report/rgn_prs_values.RData"))

quant_func <- function(layer){
  
    vec<-unlist(layer)
    q = quantile(vec,probs = c(0.025,0.25,0.5,0.75,0.975),na.rm=T)
    
    return(q)
}


#quants is a list of length 12, one for each pressure layer. Each layer has 5 quantiles
quants = map(mci_prs_vals,quant_func)
q = quants%>%
      as.data.frame()%>%
          t()%>%
          as.data.frame()%>%
          tibble::rownames_to_column()%>%
          mutate(pressure = str_replace(rowname,'_combo.tif',""))%>%
          rename(verylow = `2.5%`,
                 low     = `25%`,
                 medium  = `50%`,
                 high    = `75%`,
                 veryhigh = `97.5%`)%>%
          dplyr::select(pressure, verylow, low, medium, high, veryhigh)

df = data.frame()

for (i in 1:12){
   
     #get the pressure name
     for (j in 1:26){
       
     pres_name = names(mci_prs_vals[i])
     region   = names(mci_prs_vals[[i]])[j]
     
     #select the values
      m = mci_prs_vals[[pres_name]][[j]]
      m = m[!is.na(m)]

      h = filter(q,pressure == pres_name)
  
   n <- c(length(m[m<=h$verylow]),
          length(m[m<=h$low & m>h$verylow]),
          length(m[m<=h$high & m>h$low]),
          length(m[m<=h$veryhigh & m>h$high]),
          length(m[m>h$veryhigh]))
   
   d = data.frame(n)%>%
        mutate(perc = round(n/length(m),3),
               cat = c("verylow","low","medium","high","veryhigh"),
               brk = c(0,2.5,25,75,97.5), #assigning breaks for plotting later
               pressure = pres_name,
               rgn = region,
               mean = round(mean(m,na.rm=T),3),
               max = round(max(m,na.rm=T),3),
               total = round(sum(m,na.rm=T),3))
     
    df = rbind(df,d) 
     }
 }

Proportional area and impact of stressors across all regions

pres_names <- read_csv('pressure_names.csv')

df = df%>%left_join(pres_names,by='pressure')%>%
      mutate(rgn = ifelse(rgn == "United States Oregon", "Oregon",
                   ifelse(rgn == "United States Central California","Central California",
                   ifelse(rgn == "United States Washington", "Washington",
                   ifelse(rgn == "United States New England","New England",
                   ifelse(rgn == "United States Gulf of Mexico","Gulf of Mexico",
                   ifelse(rgn == "United States Alaska Arctic","Alaska Arctic",
                   ifelse(rgn == "Puerto Rico and Virgin Islands of the United States","Puerto Rico & US Virgin Islands",
                   ifelse(rgn == "United States Hawaii","Hawaii",
                   ifelse(rgn == "United States Northern California","Northern California",
                   ifelse(rgn == "United States Southern California","Southern California",
                   ifelse(rgn == "United States South Atlantic","South Atlantic",
                   ifelse(rgn == "United States Alaska Pacific","Alaska Pacific",
                   ifelse(rgn == "United States Mid Atlantic","Mid-Atlantic",rgn))))))))))))))%>%
      mutate(rgn_wrap = str_wrap(rgn, width = 20),
             displayName_w = str_wrap(displayName,width = 30))

#by pressure
df$cat <- factor(df$cat, levels = rev(c("verylow","low","medium","high","veryhigh")))

#climate change & pollution
cc_pol<- df%>%filter(pressure %in% c('sst','uv','ocean_acidification','inorganic','plumes_fert','plumes_pest'))
cc_pol$displayName = factor(cc_pol$displayName, levels=c('Inorganic Pollution','Organic Pollution','Nutrient Pollution','Ocean Acidification','Sea Surface Temperature','Ultraviolet Radiation'))
#fisheries
fis <- df%>%filter(pressure %in% c('artisanal_fishing','demersal_destructive_fishing','demersal_nondest_high_bycatch','demersal_nondest_low_bycatch',
                                   'pelagic_low_bycatch','pelagic_high_bycatch'))

# function to plot each of the datasets 
perc_plot <- function(data){ 
  ggplot(data) + 
    geom_bar(aes(x = rgn,y = perc,fill = factor(cat)),position = 'stack',stat = 'identity', alpha = 1) +
  scale_fill_manual(values = c('verylow' = '#4dac26',
                               'low' = '#b8e186', 
                               'medium' = '#f1b6da',
                               'high' = '#d01c8b',
                               'veryhigh'='#67001f'),
                     name = "Category",
                     labels = c("verylow","low","medium","high","veryhigh",
                    guide = guide_legend(reverse = TRUE))) +
    coord_flip()+
  labs(y = "Percent of region's area")+
  facet_wrap(~displayName_w)+
  theme(text = element_text(size = 14),
        axis.title.y=element_blank(),
        axis.text=element_text(size=14))
}

perc_plot(cc_pol)

perc_plot(fis)

#create small dataframe that assigns each pressure a short name for display purposes
s <- data.frame(pressure = unique(df$pressure),
                short    = c("art","ddf","dnd-hb","dnd-lb","inorganic","oa","pel-hb","pel-lb","fert","pest","sst","uv"))

df_dt = df%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
        dplyr::select(perc,cat,rgn,shortname)

DT::datatable(df_dt)

avg <- df%>%
       dplyr::select(pressure,rgn,rgn_wrap,mean,displayName)%>%
        unique()

#expand color palette
colourCount = length(unique(avg$pressure))
getPalette = colorRampPalette(brewer.pal(9, "Set1"))

ggplot(avg) + geom_bar(aes(y = mean, x = reorder(rgn_wrap,-mean),fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "All Pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 14),
        axis.text=element_text(size=14),
        legend.justification = c(1, 1), 
        legend.position = c(1, 1))

#without climate stressors
climate  = c("sst","uv","ocean_acidification")
avg_noclim <- avg%>%
                filter(!pressure%in% climate)

ggplot(avg_noclim) + geom_bar(aes(y = mean, x = reorder(rgn_wrap,-mean), fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "Without climate pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 14),
        axis.text=element_text(size=14),
        legend.justification = c(1, 1), 
        legend.position = c(1, 1))

avg_dt = avg%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
        dplyr::select(-pressure)%>%
        spread(shortname,mean)

DT::datatable(avg_dt)

---

British Columbia

library(tmap)
bc_map <- readOGR('shapefiles','ohibc_rgn', verbose = F)

tm_shape(bc_map)+
  tm_polygons(col = 'rgn_name',
              title = 'Regions',
              palette = 'Set2')+
tm_legend(text.size=1.5,
    title.size=1.5)

get_bc_vals <- function(layer){
  
    if(basename(layer) %in% c("artisanal_fishing_combo.tif","inorganic_combo.tif","plumes_fert_combo.tif","plumes_pest_combo.tif")){
      
       vals = raster(layer)%>%
              crop(three_ras_bc)%>%
              mask(three_ras_bc)%>%
              raster::extract(bc,method = 'simple',progress='text')%>%
              setNames(bc@data$rgn_nam)
      
    }else{
      
        vals = raster(layer)%>%
              crop(bc)%>%
              raster::extract(bc,method = 'simple',progress='text')%>%
              setNames(bc@data$rgn_nam)
    }  
return(vals)
}


#check to see if file exists.
if(!file.exists(file.path(dir_M,"git-annex/mci-report/bc_prs_values.RData"))){

  #list of all values within each region for each pressure
bc_prs_vals <- mclapply(layers,get_bc_vals,mc.cores=12)%>%
               setNames(basename(layers))%>%
               setNames(str_replace(names(.),"_combo.tif",""))

  save(bc_prs_vals, file=file.path(dir_M,"git-annex/mci-report/bc_prs_values.RData"))
}
  #vals is the list, let's try to use map to calculate mean, max
  
  load(file.path(dir_M,"git-annex/mci-report/bc_prs_values.RData"))

#quants is a list of length 12, one for each pressure layer. Each layer has 5 quantiles
bc_quants = map(bc_prs_vals,quant_func)
q = bc_quants%>%
      as.data.frame()%>%
          t()%>%
          as.data.frame()%>%
          tibble::rownames_to_column()%>%
          mutate(pressure = str_replace(rowname,'_combo.tif',""))%>%
          rename(verylow = `2.5%`,
                 low     = `25%`,
                 medium  = `50%`,
                 high    = `75%`,
                 veryhigh = `97.5%`)%>%
          dplyr::select(pressure, verylow, low, medium, high, veryhigh)

Proportional area and impact of stressors in British Columbia

df = data.frame()

for (i in 1:12){
   
     #get the pressure name
     for (j in 1:7){
       
     pres_name = names(bc_prs_vals[i])
     region   = names(bc_prs_vals[[i]])[j]
     
     #select the values
      m = bc_prs_vals[[pres_name]][[j]]
      m = m[!is.na(m)]

      h = filter(q,pressure == pres_name)
  
   n <- c(length(m[m<=h$verylow]),
          length(m[m<=h$low & m>h$verylow]),
          length(m[m<=h$high & m>h$low]),
          length(m[m<=h$veryhigh & m>h$high]),
          length(m[m>h$veryhigh]))
   
   d = data.frame(n)%>%
        mutate(perc = round(n/length(m),3),
               cat = c("verylow","low","medium","high","veryhigh"),
               brk = c(0,2.5,25,75,97.5), #assigning breaks for plotting later
               pressure = pres_name,
               rgn = region,
               mean = round(mean(m,na.rm=T),3),
               max = round(max(m,na.rm=T),3),
               total = round(sum(m,na.rm=T),3))
     
    df = rbind(df,d) 
     }
 }

df = df%>%
      left_join(pres_names,by='pressure')%>%
      mutate(rgn_wrap = str_wrap(rgn, width = 20),
             displayName_w = str_wrap(displayName,width = 30))

#by pressure
df$cat <- factor(df$cat, levels = rev(c("verylow","low","medium","high","veryhigh")))

#climate change & pollution
cc_pol<- df%>%filter(pressure %in% c('sst','uv','ocean_acidification','inorganic','plumes_fert','plumes_pest'))
cc_pol$displayName = factor(cc_pol$displayName, levels=c('Inorganic Pollution','Organic Pollution','Nutrient Pollution','Ocean Acidification','Sea Surface Temperature','Ultraviolet Radiation'))
#fisheries
fis <- df%>%filter(pressure %in% c('artisanal_fishing','demersal_destructive_fishing','demersal_nondest_high_bycatch','demersal_nondest_low_bycatch',
                                   'pelagic_low_bycatch','pelagic_high_bycatch'))

# function to plot each of the datasets 
perc_plot <- function(data){ 
  ggplot(data) + 
    geom_bar(aes(x = rgn,y = perc,fill = factor(cat)),position = 'stack',stat = 'identity', alpha = 1) +
  scale_fill_manual(values = c('verylow' = '#4dac26',
                               'low' = '#b8e186', 
                               'medium' = '#f1b6da',
                               'high' = '#d01c8b',
                               'veryhigh'='#67001f'),
                     name = "Category",
                     labels = c("verylow","low","medium","high","veryhigh",
                    guide = guide_legend(reverse = TRUE))) +
    coord_flip()+
  labs(y = "Percent of region's area")+
  facet_wrap(~displayName_w)+
  theme(text = element_text(size = 14),
        axis.title.y=element_blank(),
        axis.text=element_text(size=14))
}

perc_plot(cc_pol)

perc_plot(fis)

df_dt = df%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
        dplyr::select(perc,cat,rgn,shortname)

DT::datatable(df_dt)

avg <- df%>%
       dplyr::select(pressure,rgn,mean,displayName)%>%
        unique()

ggplot(avg) + geom_bar(aes(y = mean, x = reorder(rgn,-mean),fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "All pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 12),
        axis.text=element_text(size=12),
        legend.text=element_text(size=10))

#without climate pressures

avg_noclim <- avg%>%
                filter(!pressure%in% climate)

ggplot(avg_noclim) + geom_bar(aes(y = mean, x = reorder(rgn,-mean),fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "Without climate pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 12),
        axis.text=element_text(size=12),
        legend.text=element_text(size=10))

#create small dataframe that assigns each pressure a short name for display purposes
s <- data.frame(pressure = unique(df$pressure),
                short    = c("art","ddf","dnd-hb","dnd-lb","inorganic","oa","pel-hb","pel-lb","fert","pest","sst","uv"))

avg <- df%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
       dplyr::select(shortname,rgn,mean)%>%
        unique()%>%
        spread(shortname,mean)

DT::datatable(avg)

---

Northeast

ne_map <- readOGR('shapefiles','ohine_rgn', verbose = F)

tm_shape(ne_map)+
  tm_polygons(col = 'rgn_name',
              title = 'Regions',
              palette = 'Set2')+
tm_legend(text.size=1.5,
    title.size=1.5,
    position = c("left","top"))

get_ne_vals <- function(layer){
  
    if(basename(layer) %in% c("artisanal_fishing_combo.tif","inorganic_combo.tif","plumes_fert_combo.tif","plumes_pest_combo.tif")){
      
       vals = raster(layer)%>%
              crop(three_ras_ne)%>%
              mask(three_ras_ne)%>%
              raster::extract(ne,method = 'simple',progress='text')%>%
              setNames(ne@data$rgn_nam)
      
    }else{
      
        vals = raster(layer)%>%
              crop(ne)%>%
              raster::extract(ne,method = 'simple',progress='text')%>%
              setNames(ne@data$rgn_nam)
    }  
return(vals)
}


#check to see if file exists.
if(!file.exists(file.path(dir_M,"git-annex/mci-report/ne_prs_values.RData"))){

  #list of all values within each region for each pressure
ne_prs_vals <- mclapply(layers,get_ne_vals,mc.cores=12)%>%
               setNames(basename(layers))%>%
               setNames(str_replace(names(.),"_combo.tif",""))

  save(ne_prs_vals, file=file.path(dir_M,"git-annex/mci-report/ne_prs_values.RData"))
}
  #vals is the list, let's try to use map to calculate mean, max
  
  load(file.path(dir_M,"git-annex/mci-report/ne_prs_values.RData"))

#quants is a list of length 12, one for each pressure layer. Each layer has 5 quantiles
ne_quants = map(ne_prs_vals,quant_func)
q = ne_quants%>%
      as.data.frame()%>%
          t()%>%
          as.data.frame()%>%
          tibble::rownames_to_column()%>%
          mutate(pressure = str_replace(rowname,'_combo.tif',""))%>%
          rename(verylow = `2.5%`,
                 low     = `25%`,
                 medium  = `50%`,
                 high    = `75%`,
                 veryhigh = `97.5%`)%>%
          dplyr::select(pressure, verylow, low, medium, high, veryhigh)

Proportional area and impact of stressors in the Northeast USA

df = data.frame()

for (i in 1:12){
   
     #get the pressure name
     for (j in 1:9){
       
     pres_name = names(ne_prs_vals[i])
     region   = names(ne_prs_vals[[i]])[j]
     
     #select the values
      m = ne_prs_vals[[pres_name]][[j]]
      m = m[!is.na(m)]

      h = filter(q,pressure == pres_name)
  
   n <- c(length(m[m<=h$verylow]),
          length(m[m<=h$low & m>h$verylow]),
          length(m[m<=h$high & m>h$low]),
          length(m[m<=h$veryhigh & m>h$high]),
          length(m[m>h$veryhigh]))
   
   d = data.frame(n)%>%
        mutate(perc = round(n/length(m),3),
               cat = c("verylow","low","medium","high","veryhigh"),
               brk = c(0,2.5,25,75,97.5), #assigning breaks for plotting later
               pressure = pres_name,
               rgn = region,
               mean = round(mean(m,na.rm=T),3),
               max = round(max(m,na.rm=T),3),
               total = round(sum(m,na.rm=T),3))
     
    df = rbind(df,d) 
     }
 }

df = df%>%
      left_join(pres_names,by='pressure')%>%
      mutate(rgn_wrap = str_wrap(rgn, width = 20),
             displayName_w = str_wrap(displayName,width = 30))

#by pressure
df$cat <- factor(df$cat, levels = rev(c("verylow","low","medium","high","veryhigh")))

#climate change & pollution
cc_pol<- df%>%filter(pressure %in% c('sst','uv','ocean_acidification','inorganic','plumes_fert','plumes_pest'))
cc_pol$displayName = factor(cc_pol$displayName, levels=c('Inorganic Pollution','Organic Pollution','Nutrient Pollution','Ocean Acidification','Sea Surface Temperature','Ultraviolet Radiation'))
#fisheries
fis <- df%>%filter(pressure %in% c('artisanal_fishing','demersal_destructive_fishing','demersal_nondest_high_bycatch','demersal_nondest_low_bycatch',
                                   'pelagic_low_bycatch','pelagic_high_bycatch'))

# function to plot each of the datasets 
perc_plot <- function(data){ 
  ggplot(data) + 
    geom_bar(aes(x = rgn,y = perc,fill = factor(cat)),position = 'stack',stat = 'identity', alpha = 1) +
  scale_fill_manual(values = c('verylow' = '#4dac26',
                               'low' = '#b8e186', 
                               'medium' = '#f1b6da',
                               'high' = '#d01c8b',
                               'veryhigh'='#67001f'),
                     name = "Category",
                     labels = c("verylow","low","medium","high","veryhigh",
                    guide = guide_legend(reverse = TRUE))) +
    coord_flip()+
  labs(y = "Percent of region's area")+
  facet_wrap(~displayName_w)+
  theme(text = element_text(size = 14),
        axis.title.y=element_blank(),
        axis.text=element_text(size=14))
}

perc_plot(cc_pol)

perc_plot(fis)

df_dt = df%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
        dplyr::select(perc,cat,rgn,shortname)

DT::datatable(df_dt)

avg <- df%>%
       dplyr::select(pressure,rgn,mean,displayName)%>%
        unique()

ggplot(avg) + geom_bar(aes(y = mean, x = reorder(rgn,-mean),fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "All pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 12),
        axis.text=element_text(size=12),
        legend.justification = c(1, 1), 
        legend.position = c(1, 1),
        legend.text=element_text(size=8))

#without climate pressures

avg_noclim <- avg%>%
                filter(!pressure%in% climate)

ggplot(avg_noclim) + geom_bar(aes(y = mean, x = reorder(rgn,-mean),fill = factor(displayName)), stat = "identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(x = "Region",
       y = "Mean Impact Score",
       title = "Without climate pressures") + 
  scale_fill_manual("Pressure",values = getPalette(colourCount))+
  theme(text = element_text(size = 12),
        axis.text=element_text(size=12),
        legend.justification = c(1, 1), 
        legend.position = c(1, 1),
        legend.text=element_text(size=8))

#create small dataframe that assigns each pressure a short name for display purposes
s <- data.frame(pressure = unique(df$pressure),
                short    = c("art","ddf","dnd-hb","dnd-lb","inorganic","oa","pel-hb","pel-lb","fert","pest","sst","uv"))

avg <- df%>%
        mutate(shortname = s$short[match(pressure,s$pressure)])%>%
       dplyr::select(shortname,rgn,mean)%>%
        unique()%>%
        spread(shortname,mean)

DT::datatable(avg)